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ABSTRACT 

A two-step hybrid perturbation-Galerkin method to solve a variety of 
applied mathematics problems which involve a small parameter is presented. 
The method consists of: (1) the use of a regular or singular perturbation 
method to determine the asymptotic expansion of the solution in terms of the 
small parameter; (2) construction of an approximate solution in the form of a 
sum of the perturbation coefficient functions multiplied by (unknown) ampli- 
tudes (gauge functions); and (3) the use of the classical Bubnov-Galerkin 
method to determine these amplitudes. This hybrid method has the potential of 
overcoming some of the drawbacks of the perturbation method and the Bubnov- 
Galerkin method when they are applied by themselves, while combining some of 
the good features of both. The proposed method is applied to some singular 
perturbation problems in slender body theory. The results obtained from the 
hybrid method are compared with approximate solutions obtained by other 
methods, and the degree of applicability of the hybrid method to broader prob- 
lem areas is discussed. 
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1. INTRODUCTION 


A two-step hybrid analysis technique, which combines perturbation tech- 
niques with the Galerkin method, promises to be useful in the analysis of a 
wide variety of differential equations type problems. The hybrid technique 
was apparently first studied by Ahmed K. Noor and collaborators in conjunction 
with the finite element analysis of geometrically nonlinear problems in 
structural mechanics [22-39]. 

The Galerkin method [7] has, of course, been known and used for a long 
time. But a principal problem associated with its successful application lies 
in the choice of appropriate basis functions. It was proposed by Nagy [19] 
and developed by others [1,2,20] that the expense of analysis of certain non- 
linear structures problems can be sharply reduced by selecting a relatively 
small number of linear buckling modes as trial functions for use in a Galerkin 
technique. Then the set of solutions corresponding to a continuously varying 
parameter, such as a load factor, can be approximately determined, at least 
over some limited range of loading, with a relatively small number of non- 
linear discretized equations, in place of the much larger number of nonlinear 
finite element equations needed for modeling even relatively simple geometric 
configurations. 

In a series of papers [22-39], Noor and his collaborators have shown for 
a variety of structural mechanics problems that the terms in a Taylor series 
expansion of the solution of a parameterized system of discretized equations 
can be particularly effective as Galerkin trial functions (or basis vec- 
tors). In addition, it has been repeatedly demonstrated that the "reduced- 
basis" solutions can be useful for significantly larger values of the expan- 
sion parameter than the truncated Taylor series solutions on which they are 
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based. Noor and collaborators [29,24] have also applied the same general 
principles but without discretization to some thermal analysis and structures 
problems, A treatment of the reduced basis method from a mathematical point 
of view is given by Fink and Rheinbolt [6], 

Some general observations about the technique are the following. First, 
in many perturbation problems, much effort has to be expended to compute each 
additional term (analytically) in a perturbation expansion. Through the use 
of the proposed hybrid method, the known perturbation terms can be exploited 
more fully. Secondly, another way of viewing the technique is to recognize 
that in many perturbation expansions the functional form of the higher order 
terras can be well approximated by a linear combination of the lower order 
terms. Thus, much of the effect of the higher terras may be included by apply- 
ing the reduced basis technique to a small number of lower order terms. 
Finally, preliminary unpublished investigations indicate that, while the use 
of a Taylor series expansion is frequently limited by a finite radius of con- 
vergence, the proposed hybrid method can sometimes yield good results even 
well outside the radius of convergence. 

It is our belief that the junction of perturbation and Galerkin tech- 
niques can be useful in a wide variety of application areas, and in these 
applications the hybrid technique will give better approximations than the 
perturbation method alone. In this paper, we will present applications of the 
technique, independent of finite element or finite difference methods, and in 
an area well apart from structural mechanics. In particular, in the next sec- 
tion, we shall describe the method in more detail and then apply it to a 
simple two point boundary value problem in Section 3. In Sections 4 and 5, we 
shall apply the method to two singular perturbation problems in slender body 
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theory. In Section 6, we make some observations about our method and indicate 
some areas for further study. 


2. DESCRIPTION OF THE METHOD 

The method wish to describe is a two-step hybrid analysis technique. 
It is based upon the successive use of a straightforward perturbation expan- 
sion method and the classical Bubnov-Galerkin approximation technique. In the 
perturbation method, the solution to a particular problem involving a small 
parameter is developed in terms of a series of unknown "perturbation" func- 
tions with preassigned coefficients, i.e., gauge functions. The perturbation 
functions are usually determined by solving a recursive set of differential 
equations which are, in general, simpler than the original governing differen- 
tial equation. By contrast, in the Bubnov-Galerkin technique one seeks an 
approximate solution to the problem in the form of a linear combination of 
specified (known) coordinate functions with unknown coefficients. The coeffi- 
cients are determined by demanding that the residual formed by substituting 
the trial solution into the governing differential equations is orthogonal to 
each of the coordinate functions. 

While the perturbation and Galerkin methods are useful and have been 
successful in providing approximate solutions to a wide variety of nonlinear 
(and otherwise difficult) problems, each has certain drawbacks. The perturba- 
tion method has at least two major drawbacks. First, as the number of terms 
in the perturbation expansion increases, the mathematical complexity of the 
equations which determine the unknown functions increases rapidly. Thus, in 
most practical applications, the perturbation series is limited to only a few 
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terms. A second drawback to the perturbation method is the requirement of 
restricting the perturbation parameter to small values in order to obtain 
solutions of acceptable accuracy. (These drawbacks of the perturbation method 
have been recognized and several modifications or extensions have been 
proposed, see e.g. , Van Dyke [42] and Andersen and Geer [3].) The main short- 
coming of the Bubnov-Galerkin method is the difficulty, from a practical point 
of view, of selecting good coordinate functions. 

To illustrate the general ideas of the hybrid (or "reduced basis") 
method, suppose we are seeking (an approximation to) the solution U to the 
problem 

(2.1) L(U,e) = 0 

where L is some differential (or integral) operator and e is a small 
parameter. Here (2.1) holds in some domain D, and, in addition, U must 
satisfy certain conditions on the boundary of D. Now the application of the 
hybrid perturbation/Galerkin method can be divided into the following two dis- 
tinct steps: (1) generation of the coordinate functions by a straightforward 

perturbation expansion of U and (2) computation of the amplitudes of these 
coordinate functions by using the Bubnov-Galerkin method. 

To describe this idea in more detail, suppose that the solution to (2.1) 
can be expanded by a perturbation technique into a series of the form 

n-1 

(2.2) U = l u y .(e) + 0( Y (c)) 

j=0 J J n 

where each Uj is independent of e and {y j Ce )> is an appropriate 
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asymptotic sequence of gauge functions. The equations to determine the Uj 
are obtained by substituting (2.2) into (2.1) and setting the coefficient of 
Yj equal to zero, for j = 0,l,»»»,n-l. In a similar manner, the boundary 

conditions for each uj are determined by using (2.2) in the boundary condi- 
tion for U. 

The perturbation functions u^ are now chosen as coordinate functions 
for the Bubnov-Galerkin technique; and an approximation T) for U is 
sought in the form 


(2.3) 


n-1 

U = l u.S . 

j=0 J ^ 


where the (unknown) parameters 6^ = 6 j (e ) represent the amplitudes of the 
coordinate functions Uj. To determine these parameters, we apply the Bubnov- 
Galerkin technique to the governing equation (2.1). Thus, we substitute (2.3) 
into (2.1) and demand that the residual be orthogonal to the n coordinate 
functions over the domain D, i.e.. 


n-1 

(2.4) | L( ^ u.6.,e)u.dx = 0, k = 0, 1,» •• ,n-l. 

D j=0 J 3 K 

Equations (2,4) represent a set of n equations for the n unknown ampli- 
tudes, While (2.4) must, in general, be solved numerically, solving it is 
much simpler than numerically solving (2.1). In particular, for a fixed value 
of e, the solution to (2.4) is a point in n-diraensional space, where n is 
reasonably small, while the solution of (2.1) is a continuous function. 

We should note that this particular choice of coordinate functions over- 
comes the main drawback of the Bubnov-Galerkin method. By the way they are 
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constructed, the perturbation coordinate functions are (under certain assump- 
tions) elements of a set of functions which span the space of solutions in a 
neighborhood of their point of generation. Thus, they should fully character- 
ize the solution U in that neighborhood. Also, in many applications, the 
functions uj are determined by solving a set of linear equations, even 
though the original operator L may be nonlinear. The first property is 
necessary for the convergence of the Bubnov-Galerkin method, while the second 
property enhances the effectiveness of the proposed hybrid method for solving 
nonlinear problems. 

Another important property of the proposed method is that the coordinate 
functions, i.e., the perturbation functions, do not need to come from a 
regular perturbation expansion . In fact, all that is needed is a formal 
asymptotic expansion of the solution to (2.1) for small values of e in the 
form of (2.2), where the (Yj(e)} are a set of appropriate gauge functions, 
e.g. , expressions which involve log(e) or fractional powers of e. 
Thus , the proposed method has the potential of being applied to singular as 
well as regular perturbation problems . In the following section, we shall 
illustrate our method using the regular perturbation expansion of a simple two 
point boundary value problem, while in Sections 4 and 5 we shall apply the 
method to some singular perturbation problems in slender body theory. 


3. A SIMPLE EXAMPLE 

To illustrate the ideas just discussed, we consider the following simple 
example. Consider the two-point boundary value problem: 
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(3.1) 


U" - e U' = e , 0 < x < 1 , 


with U(0) =0 and U(l) = 0. 


The exact solution to this problem is U = (1 - e ex )/(l - e e ) - x, which 
exhibits boundary layer behavior around x = 1 for large values of e. A 
regular perturbation expansion of U yields a series of the form 


(3.2) 


U = l u.(x)e J + 0(e n+ *) 
j=l J 


where each Uj is a polynomial in x of degree (j+1) and vanishes at x « 
0 and 1. Following the ideas discussed above, we look for an approximate 
solution U in the form 


(3.3) 


U =■ l u. (x) 6 .(e) 
j=l 3 3 


where the amplitudes 6^ are to be determined. To determine them, we sub- 
stitute (3.3) into (3.1) and apply the Bubnov-Galerkin criterion (2.4) to ob- 
tain the set of equations 


(3.4) 


I 3j ^(e ) 6^ — bj(c), j — 1 ,2,» • • ,n 


where 


• = J u.(x)[u~(x) - eu'(x)]dx, b =e / u (x)dx. 

J 9 K * J n ■ J 
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Ia this case, the equations for the 6^ are linear since the original 

problem is linear. For example, setting n = 2 in (3.4), we find 


1 l+e 2 /60 


, s 2 


1+e /60 


For small values of e , we see that these expressions reduce to e and 
o 

e , as they should. For n = 4, we find (by using the symbolic computation 
system MACSYMA [16]) that 


. _ 15120e + 420e 

6 1 2 ZT 

15120 + 420e + e 


, 6 2 = e6j 


*3 " 


15120e 

2 2T ’ 

15120 + 420e + e 


64 eSy 


(Note that the S^'s remain finite as e becomes large, which suggests 
that the approximate hybrid solution for this problem may be valid even for 
large values of e . ) 

In Figures 1 and 2, we illustrate the comparison of the exact solution 
with our approximate solution and the perturbation solution for n = 2 and n 
= 4, respectively. For n « 2, the three solutions agree very well for values 
of e up to about e » 1. For larger values of e, the perturbation 
solution begins to diverge from the exact solution, while U still gives 
reasonable approximations up to about e = 5. A similar comparison using 
four terras in the approximate solutions shows that our hybrid solution is con- 
sistently more accurate than the perturbation solution and gives reasonable 
results, even within the boundary layer, for values of e up to about 10, 
where the perturbation solution is meaningless. 
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Sorae of the reasons why our hybrid method provides a good approximation 
in this case, as well as insights this example provides for the validity of 
our method, will be discussed in Section 6. The application of our hybrid 
method to several general classes of two point boundary value problems, as 
well as various extensions of the method, are currently under investigation 
and will be reported elsewhere. 

4. POTENTIAL FLOW PAST A SLENDER BODY OF REVOLUTION 

Slender body theory has been developed and applied largely in the context 
of classical fluid mechanics. For example, potential flow past an axially 
symmetric slender body has been studied extensively since the work of Munk 
[18], This work is now a part of slender body theory, which is the theory of 
fields about a slender body. It is discussed in detail, along with relevant 
references, in various books, e.g., Thwaites [40] and Van Dyke [41], as well 
as the paper by Newman [21]. More recent work on the axially symmetric flow 
problem has been done by Moran [17], Landweber [15], and Handelsman and Keller 
[12]. Recently, these results have been extended and generalized by Geer 
[8], In addition, the ideas of uniform slender body theory have been applied 
successfully to problems in other areas, such as Stokes flow past a slender 
body (Geer [9]), the electrostatic field about a dielectric slender body 
(Barshinger and Geer [4]), the scattering of an arbitrary acoustical wave by a 
slender body of revolution (Geer [10]), as well as a special case of electro- 
magnetic scattering (Geer [11]). 

To illustrate the applicability of our hybrid method to problems in this 
area, we shall consider two problems which involve a slender body of revolu— 
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tion. These problems will serve as model problems for the application of our 
method to more general problems in slender body theory. Typically, these 
problems are not regular perturbation problems, since the solutions contain 
terms involving log(e) as well as powers of e, where e is the 
slenderness ratio of the body (see e.g. , Geer [8]). In addition, the solu- 
tions to some of these problems involve series in inverse powers of log(e), 
which converge (at best) very slowly, even for reasonably small values of e. 
The specific problems we shall consider are: (1) potential flow past a 

slender body, which will be considered in this section and (2) the electro- 
static potential around a perfectly conducting slender body (Section 5). Each 
of these problems contains features representative of larger classes of prob- 
lems for which we believe the hybrid method will be of some significant 
value. For example, the flow problem is a good representative of a singular 
(as opposed to regular) perturbation problem involving a relatively rapidly 
converging series. Similar series occur in the study of certain Stokes flow 
problems (Geer [9]), acoustical scattering by a hard body (Geer [10]) or a 
perfectly conducting body (Geer [11]), the potential field about a dielectric 
body (Barshinger and Geer [4]), and certain problems involving oblate bodies 
(Horaentcovschi [14]). The electrostatic problem involves a sequence of very 
slowly converging series. Similar slowly converging series also occur in 
investigations in several other problem areas, including certain other Stokes 
flow problems (Geer [9]), scattering by a soft body (Geer [10]), and problems 
involving the interaction of ships in shallow water (Davis and Geer [5]). For 
simplicity, we shall consider here only a simple example in each of these 
areas. These model problems will serve to illustrate the usefulness of our 
method, and it will also be clear how our method can be applied to more 
general problems in these areas. 
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In this section, we shall consider the problem of determining the 
potential flow around a slender body of revolution which is immersed in a 
uniform stream parallel to its axis (see e.g. , Handelsman and Keller [12] or 
Geer [8], where this problem is generalized). We let the body be described 
by r = (x 2 + _ e ^s( z ) , for 0 _< z _< 1 (see Figure 3). Here e 

is the slenderness ratio of the body, defined as the ratio of the maximum 
radius of the body to its length, while S(z) describes the particular shape 
of the body. We assume that S(0) = 0 = S(l), while S'(0) ^ 0 ^ S'(l), where 
the primes denote differentiation with respect to z. These last conditions 
ensure that the body has a rounded (i.e., not sharp) nose and tail. (See Geer 
[8] for a discussion of these conditions.) Then, following the ideas of 
Handelsman and Keller [12], we represent the part of the velocity potential $ 
due to the presence of the body as the superposition of potentials due to 
point sources distributed along a portion of the axis of the body and lying 
entirely within the body. Thus, we write: 

g 

(4.1) #(r 2 ,z) - z - I 1 — f(g,e)dg, 

“ /(z-5) 2 +r 2 

where 4irf is the unknown density of the source distribution and a 
and 8, which determine the extent of the distribution within the body, 
satisfy 0 < a < 8 < 1 (see Figure 3). The expression (4.1) satisfies 
Laplace's equation outside the body and reduces to z at infinity. The con- 
dition that the normal component of the velocity vanish on the surface of the 
body, when used with (4.1) becomes 

g 

ye S'(z) = (^-) J " z _ ' f(C,e)dg, 

• /(z-t)WsM 


(4.2) 
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which must hold for 0 < z < 1, Equation (A. 2) is an integral equation for 
the source density f(z,e), from which a and g can also be deter- 

mined. In particular, Handelsman and Keller have shown that f can be deter- 
mined as a series of the form 

00 i- 1 

(4.3) f (z ,e ) = l e 2i l (log( e 2 )) j f. ,(z) 

i=l j=0 ’ 3 

where the coefficient functions f^ j(z) can be determined recursively. In 
particular, using some of the results of Handelsman and Keller, we can write 

2 

f(z,e) = S'(z) + e 4 log(e 2 )f 2 j ( z) + e 4 f 2 Q (z) 

(4.4) + O(e 6 (log(e 2 )) 2 ) , 

£ 2>! ( 2 > - ^ (s(os f 2 j 0 (z) - £ 

where is a certain linear operator, defined precisely by Handelsman 

and Keller (see also Geer [9], appendix B). The parameters a and S 

2 

can be determined as power series in e and are given by 

a = C^) 2 - C 1 C 2 (|) 4 + + 2C 2 )(|) 6 + 0(e 8 ), C. = S <j) (0)/j! 

(4.5) 

3 = 1- d x (|) 2 + d 1 d 2 (|) 4 - d 1 (d 1 d 3 + 2d 2 )(|) 6 + 0(e 8 ) , df (-l) j S (j) (l)/j ! . 

To illustrate the application of the hybrid method to this problem, we 
shall treat a and g as known quantities (since they can be expressed 

as Taylor series in e , with several coefficients having been determined), 
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aad concentrate on the behavior of f with e. Following the ideas pre- 
sented in Section 2, we look for an approximation 1: to f in the form 

n i-1 

(4.6) ^ ( z ,e ) = l l f. .(z) 6. .(e) 

i-0 j=0 1,J 1J 

where the coefficients f^ ^(z) are those which appear in (4.3) above and the 

amplitudes 6. .(e) must be determined. To determine them, we substitute 
1 > 3 

(4.6) into (4.2), apply the orthogonality conditions, and obtain the condi- 
tions: 


(4.7) 


n i-1 

Vl = 5 l 0 ^,1,1,3 6 i,j’ 


where 


2 1 


"m** / s ' (z > f k,i (z)d2 


nj-I %/ 


1 d (z - O 


a / (z-£ ) 2 +e 2 s(z) 


f i,j (5)d? > f k,l (z)dz - 


Equations (4.7) are a system of linear equations for the unknown ampli- 
tudes 6^, which can be solved in a straightforward manner. As a special 
example, we consider the case when n = 1 and hence our approximate solu- 
tion t is of the form 


(4.8) 


V T- jV<0)V c i TF 3F / 

’ 0 ’ * ’ 0 a 


f(z,e) - •£• S'(z)5 1>Q , 6^ 0 b i > 0 yrc l, 0,1,0’ 

i , e 




a / (z-C) 2 +c 2 s(z) 


S'(£)d£dz. 
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From (4.8) we can easily show that 


6^ g * + Y^e^log(e^) + y 2 ^ + O(e^(log(e^))^» as e ■*■ 0, 


where and Y 2 are certa i- n constants. Hence, as e 0, our 
approximate solution (4.8) reproduces the first term in the perturbation solu- 
tion (4.4) exactly, while also "anticipating" some of the e behavior of 
the higher order terras. For the special case when the body is a prolate 
spheroid, i.e. , when S(z) = 4z(l-z), the integrals appearing in (4.8) can be 
evaluated explicitly, yielding 


(4.9) 


6j Q - e 2 [/ 1 - 4 e 2 - 2 £ 2 C(e)] 


C(e ) = 21og[— 


vW 

Tz 


When (4.9) is inserted into (4.8), we find that (4.8) yields the exact solu- 
tion to the integral equation (4.2) for this case. Thus, the hybrid method, 
for this special case, yields the exact solution using only the "inf ormation" 
contained in the first term of the perturbation solution. This point will be 
discussed in more detail in Section 6. 

To obtain some idea of the accuracy of our approximate solution (4.6) for 
other body shapes, we have applied our method to a class of "dumbbell" shaped 
bodies, described by S(z) = 4bz(l-z){l-bz(l-z)}, with 2 < b < 4. In Figure 
4, we illustrate our results for the representative case b « 3 and for e = 
0.03, 0.07, and 0.11. In particular, we compare our hybrid solution (4.6) 
with the perturbation solution (4.3) and a solution to (4.2) obtained by pure- 
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ly numerical means (see the Appendix). As the figure clearly illustrates, the 
hybrid solution gives consistently better results than the perturbation solu- 
tion alone, especially for the case e = 0.11, which appears to be close to 
the radius of convergence of the perturbation solution. In addition, we have 
used our hybrid solution in (4.1) to calculate the (non-dimensional) pressure 
coefficient C^(z) = 1 “ |^$ | on the surface of a slender body. In Figure 

5, we have plotted for the prolate spheroid with e = 0.3, using our 
hybrid solution, the perturbation solution, and the exact solution. In Figure 

6, we have made a similar comparison for our dumbbell shaped body with b = 3, 
which was also discussed by Wong, Liu, and Geer [43]. Figure 5 illustrates 
the facts that our hybrid method has reproduced the exact solution to the flow 
problem and that the perturbation solution is (slowly) converging to the cor- 
rect result. Figure 6 again illustrates that the hybrid method produces con- 
sistently better results than the perturbation solution alone. 


5. ELECTROSTATIC POTENTIAL FIELD ABOUT A SLENDER CONDUCTOR 

We now consider the problem of determining the electrostatic potential 
field about a perfectly conducting slender body of revolution. For simpli- 
city, we again consider only a simple version of this classical problem, name- 
ly, the case when there is no external field and the body is raised to unit 
potential. Thus, using the notation of the previous section, we represent the 
potential as 


$(r ,z) = / 


° /(z-?) 2 +r 2 


f(z,e )dg 


(5.1) 
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where 4irf(z,e) is again the (unknown) source strength density. The ex- 
pression (5.1) satisfies Laplace's equation outside the body and vanishes at 
infinity. The boundary condition on the surface of the body leads to the con- 
dition 

(5.2) / 1 f(C,e)dC - 1, 0<z<l 

a /(z~C ) 2 +e 2 s(z) 

which is an integral equation for f(z,e), from which a and 3 can 
also be determined. In particular, Handelsman and Keller [13] have shown 
that a and 3 are again given by the expressions (4.5), while f has 
an asymptotic expansion of the form 

(5.3) f (z,e ) - l l e 2m (log( e 2 ))' j f ,(z) 

m=0 j=l “» J 

where the coefficient functions f m> j(z) can be determined recursively. In 
particular, using some of their results (see also Barshinger and Geer [4], 
where these results have been generalized to the case of a dielectric body) we 
can write 


f 0 j+1^ = f 0 j( z ) lo g{4z(l-z)/S(z)} 

1 

+ / {f .(c) - f n ,(z)}/|s - z| dc, j = 1 > 2 , • • • 

Q U >J 


(5.4) 
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f l j+1^ = f l j( z ) lo g{4z(l-z)/ s (z)} 

1 

+ I { f i ^5) “ f, .(z)}/|? ~ z| d£ 

Q A 9 J 1 9 J 

+ (l/4)S(z)f~ (z) + Lf .(z), j = 1,2,..., 

0,j+l U,J 

where L is a certain linear operator defined explicitly by Handelsman and 

Keller [13] (see also Geer [9], Appendix B). 

We note that the series (3,3) contains a sequence of (at best) slowly 

2 

converging series involving inverse powers of log(e )• We have used the 
expressions (5.4) to compute some of the f m j explicitly and then used these 
results in (5.3) to find approximations for the actual singularity density 
f(z,e). In Figure 7, we have plotted the approximation for f obtained from 
(5.3) using different numbers of terras* These approximations are compared 
with a solution to (5.2) obtained by purely numerical methods (see the Appen- 
dix). As can be seen from the figures, several terms are needed in the 
perturbation solution to obtain a good approximation, even for reasonably 
small values of e. 

We now apply our hybrid method to this problem and look for an approxi- 
mate solution f in the form 


(5.5) 


*<*.«>- 1 r £ ,<*>« ,(*■> 


m=0 j=l 




m » J 


where the coefficients f i(z) are those which appear in (5.3) and the 

J 

amplitudes 6 . are to be determined. To determine them, we substitute 

(5.5) into (5.2), apply the orthogonality conditions, and obtain a set of 
equations of the form (4.7), where 
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(5.6) 


3 k,l = / f k,l (z)dz 


i e 


■ i i j =/ / 


0 ° /(z-5) 2 +e 2 s(z) 


f i j(C)d?f k 1 (z)dz. 


Hence, the 6 , again satisfy a system of linear algebraic equations, 
which can be solved in a straightforward manner. In particular, for the 
special case when n = 0 and J 0 = 1 (i*e., f contains only one term), we 


find that 


(5.7) 



- {/ / t(z " 5 ) 2 + e 2 s(z)]" 1/2 d5dz) 

0 a 

(loge 2 )" 1 


(logs 2 ) 2 / log [ — j— — ] dz 


0((loge 2 ) 3 ) + O(e 2 (loge 2 ) 2 . 


Thus, 6 q ^ reproduces the first term in the asymptotic expansion (5.3) 
exactly, while again "anticipating" some of the higher order e behavior of 
the solution. For the special case when the body is a prolate spheroid with 
S(z) = 4z(l-z), we find that the integral appearing in (5.7) can be evaluated 
explicitly to yield 1 = ^ > where C(e ) is defined in (4.9). 
When this result is inserted into (5.5), we find that (5.5) yields the exact 
solution to the integral equation (5.2). These points will be discussed more 
fully in Section 6. 

In Figure 7, we have also plotted the approximation to f(z,e) ob- 
tained from our hybrid solution (5.5). As the figures clearly indicate, the 
hybrid method consistently produces better approximations than the perturba- 
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tion solution alone. Some of the reasons for these better approximations will 
be discussed in the next section. 


6. DISCUSSION AND OBSERVATIONS 

We can now make some observations about the hybrid method and briefly 
discuss, at least heuristically, why it seems to work as well as it does. 

If we return to equation 2.2, we note that the perturbation method seeks 

a 

to find an approximate solution U to the true solution U of (2.1) in the 
form of a finite linear combination of the perturbation coordinate functions 
ui with specified amplitudes (gauge functions) {y.(e)>. Our hybrid 
method, on the other hand, seeks an approximate solution U in the form of 
equation (2.3), which can be viewed as a linear combination of some of the 
perturbation coordinate functions with amplitudes {6^} which are not 
specified a priori. Instead, the {6^} are determined so that tl will 
be "the best" approximation of the form (2.3) using the Galerkin criterion 
(2.4). Thus, in a geometrical sense, the perturbation method seeks to approx- 
imate the shape of the true solution with a specific linear combination of 
the Uj , while the hybrid method seeks to approximate the shape of U with 
the "best possible" (in the Galerkin sense) linear combination of a specified 
subset of the {uj}. Thus, the perturbation method ignores, for example, 
linear dependences among the {uj}, while the hybrid method implicitly takes 
these dependences into account. An extreme case of this occurs in the prob- 
lems in Sections 4 and 5 in the special case when the body is a prolate 
spheroid. In this case, all of the perturbation coordinate functions are 
constant multiples of the first function and hence are all linearly dependent 
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(see, for example, Wong, Liu, and Geer [43]). The perturbation method thus 
produces an Infinite number of terms in a formal asymptotic series, but each 
term has the same "shape" as the first term, with only the amplitudes of the 
terms differing. The hybrid method in essence "seeks out" (in the Galerkin 
sense described) the total dependence of the shape of the solution on the 
terms included in the approximation. Thus, in the case of the prolate 
spheroid, since the exact "shape" of the true solution is contained in the 
first term of the approximation, the hybrid method produces the exact solu- 
tion. Only the amplitude of the term needs to be modified. 

The hybrid method, then, would seem to work well when the perturbation 
coordinate functions {uj} exhibit a significant degree of linear dependence. 
To illustrate this dependence for the examples we have presented, in Figures 
8-10 we have plotted several of the perturbation coordinate functions used in 
the hybrid approximation. In Figure 8 , we have plotted the first 15 perturba- 
tion functions for the simple example discussed in Section 3 (with each func- 
tion normalized so that its maximum absolute value is 1). We recall from 
their definition that each uj is a polynomial of degree j +1 and hence they 
are all independent. However, as can be seen clearly from the figure, on the 
interval [ 0 , 1 ] these functions are highly dependent, with essentially only a 
few different "shapes" occurring among the first 15 functions. We have ex- 
pressed each of these 15 funtions in terms of the 4 functions obtained by 
applying the Gram-Schmidt orthonormalization process to the first four func- 
tions uj to U 4 . In Table I, we have presented the L 2 “ norm the 
difference of these expressions and the original perturbation coordinate func- 
tions. The norms confirm that fact indicated by Figure 8 that essentially any 
linear combination of the first 15 coordinate functions can be well expressed 
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in terras of just the first four Uj* Thus, we can say that the higher order 
perturbation coordinate functions add only a small amount of new information 
about the shape of the true solution. Presumably, these first four functions 
contain much of the "shape information" contained in the perturbation func- 
tions of order higher than fifteen as well. (We have examined several other 
model two point boundary value problems and have found conclusions similar to 
those just expressed. These examples, as well as other applications of our 
method, are currently under investigation and will be reported elsewhere.) 
These observations about the linear dependency of the perturbation functions 
help to account for the accuracy of our hybrid solution. If, on the other 
hand, the perturbation functions in a particular problem turn out to be only 
weakly dependent on the domain of interest, then we would not expect our 
hybrid solution to be much of an improvement over the perturbation solution. 

In Figure 9, we have plotted the first four (normalized) perturbation 
functions used in the potential flow problem discussed in Section 4. In Table 
II, we have recorded certain I^-norms indicating the accuracy of representing 
these functions, as well as the "exact" solutions for different values of e, 
in terms of the first four coordinate functions. As the table clearly indi- 
cates, the fourth perturbation function adds very little "new" information 
about the shape of the solutions. In a similar manner, the first 10 perturba- 
tion functions used for the electrostatic problem discussed in Section 5 are 
presented in Figure 10, while Table III contains certain L 2 -norras which 
indicate the degree of dependency among the different functions and the 
"exact" solutions. Again the degree of dependency that does exist helps to 
explain why the hybrid method produces considerably better accuracy than the 
perturbation method alone. (We have investigated several other slender body 
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shapes in addition to the family of "peanut" shaped bodies presented here and 
have again found results concerning dependency similar to those indicated 
above. ) 

The ability of the hybrid method to pick out certain dependencies in 
higher order terms can be seen explicitly in another way. To illustrate this 
way for the problem considered in Section 5, we note that we can construct a 
hybrid approximate solution of the form (5.5) with only one term. When this 
is done, the single amplitude 6 q ^ is given by (5.7). Using (5.4) and 
(5.7), we note that the hybrid solution can be expressed as 

f(z,e) = f o,i^ z ^ 6 o,:^ e ^ 

(6.1) = f Q ^z){ (log(e 2 )) 1 + (/ f Q 2 (z) f o 1 (z)dz)(log(e 2 )) 2 

+ 0((log( e 2 ))" 3 )}. 

The first term in the second line of (6.1) is just the first term in the 
perturbation expansion (5.3). The second term, however, represents the 
component of the second perturbation function fg 2 which is parallel to 
fp^j. Similar observations hold for the other problems we have considered 
when the 6^'s are ex P an ded for small values of e. Based upon this 

discussion, there is another way of characterizing those problems for which 
the hybrid method should be helpful. If the perturbation method produces an 
approximate solution which seems to have about the right shape, but not the 
proper amplitude, then the hybrid method should be applied. This was very 
evident in the computation of the pressure coefficient for a prolate spheroid 
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cons idered in Section 4 (see Figure 5). A more interesting example was con- 
sidered by Davis and Geer [5], who used uniform slender body theory to 
investigate the interaction of ships maneuvering in shallow water. Their 
analysis produced expressions for the forces and torques acting on the ships 
which had about the right shape (when expressed as a function of the distance 
between the ships and compared with some experimental evidence) , but not the 
proper amplitude. Thus, it would seem plausible to expect that the hybrid 
method should improve the accuracy of their results for this problem. 

The convergence of the hybrid method as the number of terms in the ex- 
pression increases has been investigated experimentally, but has not been 
addressed formally. For instance, for our simple example of Section 3, we 
have computed the relative 1^ error of both the perturbation and hybrid 
solutions for values of e between 0 and 15, using either n = 2 or n = 
4 terms in the approximation. The results are plotted in Figure 11. Clear- 
ly, for each value of e , the error in the hybrid solution is diminishing 
as n increases from 2 to 4. Asymptotically, as e approaches infinity, 

these errors approach 0.9258 (n = 2) and 0.9045 (n ■ 4). It is interesting to 
note that the perturbation solution has a finite radius of convergence 
R = 2x, since the exact solution has singularities at e = ±2xi. (The 
denominator in the expressions for the deltas with n * 4 has a coots at 
e = ±6.306i, which are good estimates of these singularities.) The hybrid 
solution provides a reasonable approximation for values of e beyond the 
radius of convergence of the solution upon which it is based, and presumably 
this approximation will improve as n is increased. This point will be dis- 
cussed at greater length elsewhere, when we apply our method to several other 
two point boundary value problems. 
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There are several ways in which the basic method we have presented can be 
refined or extended. For example, a set of functions other than the set 
{uj} can be used in the orthogonality conditions (2.4) to determine the 
amplitudes {6^}. That is, we could demand that the residual formed by 
inserting (2.3) into (2.1) should be orthogonal to another set of functions, 
say {v.}, which might be related, say, to the adjoint of the operator L . 
In addition, some other criterion, in place of (2.4), might be used to deter- 
mine the {6^}, e.g., a least squares criterion. For computational pur- 
poses, there may be some advantage to orthogonalizing and scaling the {u^} 
before applying the conditions (2.4). Also, there is no reason why the 
first n terms in the perturbation solution must be used as a basis for the 
hybrid solution. In fact, it is certainly conceivable that certain higher 
order perturbation functions might add more new shape information about the 
solution than some lower order terms. In the problem considered in Section 5, 
for example, we could consider forming a hybrid solution consisting of, for 
example, f 0> i, f 0,2> f 0,3> and f l,2 # The P r °hlem of choosing an "optimal" 
set of perturbation functions to use as a basis for the hybrid solution (both 
from a computational and theoretical point of view) is currently under invest- 
igation and progress on this point, as well as those mentioned above, will be 
reported as progress is made. 
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APPENDIX 


In this appendix, we shall describe briefly the iterative method we used 
to solve equations (4.2) and (5.2) numerically. Using a result of Handelsman 
and Keller [12], we see that the right side of equation (4.2) approaches 
2f(z,t.) as c -*• 0; Thus, if we add and subtract 2f(z,e) on the right 
side of (4.2) and then define 

(A. 1 ) F(z,e) = 4f(z,e)/e 2 , 

we find that (4.2) can be written as 

(A. 2) S'(z) = F(z,e) + L(F(z,e)), 

where the linear operator L is defined by 

8 . _ t 

(A. 3) L(F(z) ) = (1/2) {(d/dz) / — F(t) dt - 2F(z)>. 

a /(zt) 2 +e 2 s(z) 

2 

From Handelsman and Keller [12], it follows that L(F(z)) = 0(e F(z)) as 
e 0. This suggests that we could solve equation (A. 2) by defining the 
sequence of approximations {F^(z)} by 


(A. 4) F Q (z) = S'(z) 


(A.5) 


F fc+1 (z) = L(F k (z)), k = 0,1,2, 


• • • 
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The iterative scheme (A.3)-(A.5) was used to solve equation (4.2) numerical- 
ly. As indicated in Section 4, the limits of integration a and g were 
treated as known quantities, using the Taylor series expansions obtained by 
Handelsman and Keller, as well as several more terms which we calculated using 
the symbolic computation system MACSYMA [16]. The terms in these series were 
used with the ratio and root tests to estimate the radius of convergence of 
the formal perturbation series, as mentioned in Section 4. We found that the 
scheme (A.3)-(A.5) converged rather quickly for a £z _< g, but converged 
rather slowly for 0 ^ z < a or g < z _< 1. Nevertheless, the scheme did 
produce (eventually) the correct solution to equation (4.2) on the entire 
interval 0 _< z _< 1 for all the examples we tried. 

Equation (5.2) can also be solved by iteration by first rewriting it as 

(A. 6) 1 = G(z,e)f(z,e) + M(f(z,e)), 

(A.7) G(z,e) - / 1 dt - 

a /(z-t) 2 +e 2 s(z) a-z+/ (a-z) 2 +e 2 s(z) 

(A. 8) M(f (z,e )) = J 6 f jh e 2 . Z , f . ( £>i2 dt. 

a /(z-t) 2 4e 2 s(z) 

Using a result from Handelsman and Keller [13], we see that 

G(z,e) = 0(log(e^)) and M(F(z)) = 0(F(z)), as e 0. These facts 

suggest that we define a sequence of approximations by 


(A. 9) 


£q(z » e ) = 1/G(z,e), 
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(A. 10) f k+1 (z,e) = {1 - M(f k (z,e))}/G(z,e), k = 0,1,2,... . 

We found that the scheme (A.9)-(A.10) converged to the solution of (5.2) for 
all the cases we tried, but with the same type of non-uniform convergence as 
noted above. 
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Table I 


n 


u(z) 

1 

2 

3 

4 

Uj (z): 





3 = 1 

0 

0 

0 

0 

2 

1.000 

0 

0 

0 

3 

0.189 

0.189 

0 

0 

4 

1.000 

0.100 

0.100 

0 

5 

0.225 

0.225 

0.014 

0.014 

6 

1.000 

0.123 

0.123 

0.005 

7 

0.233 

0.233 

0.018 

0.018 

8 

1.000 

0.129 

0.129 

0.007 

9 

0.235 

0.235 

0.020 

0.C20 

10 

1.000 

0.130 

0.130 

0.008 

11 

0.235 

0.235 

0.020 

0.020 

12 

1.000 

0.130 

0.130 

0.008 

13 

0.235 

0.235 

0.020 

0.020 

14 

1.000 

0.131 

0.131 

0.008 

15 

0.235 

0.235 

0.020 

0.020 

U (e = 2) 

0.122 

0.013 

0.001 

0.000 

U (e = 5) 

0.269 

0.067 

0.014 

0.003 

U (e = 10) 

0.403 

0.164 

0.062 

0.021 


n 

L 2 - norms of the differences u(x) - £ a.g.(x) for the simple two- 

j=l J J 

point boundary value problem of Section 3, where the functions u(x) have 
been normalized and the functions gj(x) have been obtained from the 
perturbation functions u^ to u^ by the Gram-Schmidt process. The numbers 
appearing in the table are the quantities 


(1 


n 


l 

3-1 


2 . 1/2 

a. ) 


where 


1 

a. - j u(z)g (z)dz. 

3 0 3 


Thus, the smaller the entry the better the function is approximated. The last 
three rows in the table correspond to the exact solution of the problem (3.1). 
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Table II 


n 


f (z) 

1 

2 

3 

4 

f l,0 

0 

0 

0 

0 

f 2, 1 

0.647 

0 

0 

0 

f 2*0 

0.975 

0.045 

0 

0 

f 3, 2 

0.885 

0.445 

0.095 

0 

f (c = .03) 

0.045 

0.000 

0.000 

0.000 

f (c = .07) 

0.195 

0.032 

0.000 

0.000 

f (c = .11) 

0.525 

0.276 

0.032 

0.032 


L, 2 ~norms of the differences f(z) - £ a.g.(z), where the functions 

j=l 3 J 

f(z) have been normalized and the functions gj(z) bave been obtained from 
the potential flow perturbation coefficient functions fj^Q, ^2,0’ aru * 

f 3 2 by the Gram-Schmidt process. The numbers appearing in the table are the 
quantities 

O o 1 /O ^ 

(1 “ l a . ) , where a = / f(z)g (z)dz. 

j=l J J 0 J 


n 


The last three rows in the table correspond to the solution of the integral 
equation (4.2) obtained numerically. 
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Table III 


n 


f(z) 

1 

2 

3 

4 

f 0 ,j< z > : 






0 

0 

0 

0 

2 

0.823 

0 

0 

0 

3 

0.985 

0.032 

0 

0 

4 

0.990 

0.089 

0.000 

0 

5 

0.991 

0.158 

0.000 

0.000 

6 

0.991 

0.245 

0.000 

0.000 

7 

0.992 

0.349 

0.032 

0.000 

8 

0.993 

0.468 

0.055 

0.000 

9 

0.994 

0.592 

0.084 

0.000 

10 

0.995 

0.706 

0.126 

0.000 

11 

0.997 

0.799 

0.179 

0.000 

12 

0.998 

0.868 

0.239 

0.000 

13 

0.998 

0.915 

0.302 

0.000 

14 

0.999 

0.946 

0.370 

0.045 

15 

0.999 

0.965 

0.439 

0.071 

f (e - 0.3) 

0.045 

0.000 

0.000 

0.000 

f (e = .07) 

0.195 

0.032 

0.000 

0.000 

f (e « . 11 ) 

0.525 

0.276 

0.032 

0.032 


n 

Lo-norms of the differences f(z) - £ a.g.(z), where the functions 

j=l J J 

f(z) have been normalized and the functions gj(z) have been obtained from 
the electrostatic potential perturbation functions ^o,l» ^ 0 , 2 * ^0 3 » an< * 

fQ 4 by the Gram-Schmidt process. The numbers appearing in the table are the 
quantities 

n 2 1 

(1 - £ a.), where a. = / f(z)g.(z) dz. 

j=l J J 0 J 


The last three rows in the table correspond to the solution of the integral 
equation (5*2) obtained numerically. 
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FIGURE CAPTIONS 


1. Comparison of solutions to the two point boundary value problem (3.1): 

(o o o o o) exact solution, (- - ) perturbation solution, ( ) 

hybrid solution, using two terms in the approximate solutions for (a) 
e = 2 and (b) e = 5. 


2. Comparison of solutions to the two point boundary value problem (3.1): 

(o o o o o) exact solution, (- - - — ) perturbation solution, ( ) 

hybrid solution, using four terras in the approximate solutions for (a) 
e = 5 and (b) e = 10. 


3. A typical slender body of revolution, with an indication of the 
coordinate system and the region of distribution of singularities. 

4. Comparison of approximate solutions to the integral equation (4.2) 

obtained numerically (o o o o o), using the perturbation solution (- 

- -) , and using the hybrid solution ( ) for: (a) e = 0.03; (b) 

e = 0.07; (c) e = 0.11. In (a), only one term was used in the 

perturbation and hybrid solutions, while in (b) and (c) three terras were 
used. 


5. Pressure coefficient for a prolate spheroid with e = 0.3 using the 
exact solution (o o o o o), the perturbation solution (-----) with 

one, three, and six terms, and the hybrid solution ( ) with one 

term. 


6. Pressure coefficient for a peanut shaped body of section 4 using the 
solution obtained numerically (o o o o o), the perturbation solution (- 
----), and the hybrid solution ( ) for (a) e = 0.03, (b) 

e ■ 0.07, and (c) e = 0.11. In (a), only one term was used in the 
perturbation and hybrid solutions, while in (b) and (c) three terms were 
used. 


7. Comparison of approximate solutions to the integral equation (5.2) 

obtained numerically (o o o o o), using the perturbation solution ( 

- -), and using the hybrid solution ( ) for: (a) e = 0.03; (b) 

e = 0.07; and (c) e = 0.11. In each case, four terms were used in 
the perturbation and hybrid solutions. 

8. The first 15 perturbation coordinate functions u^(x) appearing in 

equation (3.2). The functions u^ to u^ are plotted with solid 

lines, while the remaining functions appear as dotted lines. Each 

function has been normalized so that its maximum absolute value is 1. 
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9. The first four potential flow perturbation coordinate functions 
f^j(z) appearing in equation (4.3). The functions have been 

normalized so that each one is equal to 1 at z = 0. 

10. The first 10 electrostatic potential perturbation coordinate functions 
fQ j(z) appearing in equation (5.3). The functions have been 

normalized so that each one is equal to 1 at z = 0. 


11. Plot of the relative L 2 error between the exact solution to problem 
(3.1) and the perturbation solution (- - - - ~) and the hybrid solution 
( ) using 2 and 4 terms in these approximate solutions. 
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1(a) 
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1(b) 
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2(a) 
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Figure 7(a) 
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